Ab initio study of the atomic motion in liquid metal surfaces: comparison with 

Lennard-Jones systems. 
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It is established that liquid metals exhibit surface layering at the liquid-vapor interface, while 
dielectric simple systems, like those interacting through Lennard-Jones potentials, show a monotonic 
decay from the liquid density to that of the vapor. First principles molecular dynamics simulations 
of the liquid- vapor interface of several liquid metals (Li, Na, K, Rb, Cs, Mg, Ba, Al, Tl and Si), 
and the NaaKr alloy near their triple points have been performed in order to study the atomic 
motion at the interface, mainly at the outer layer. Comparison with results of classical molecular 
dynamics simulations of a Lennard-Jones system shows interesting differences and similarities. The 
probability distribution function of the time of residence in a layer shows a peak at very short 
times and a long lasting tail. The mean residence time in a layer increases when approaching the 
interfacial region, slightly in the Lennard-Jones system, but strongly in the metallic systems. The 
motion within the layers, parallel to the interface, can be described as diffusion enhanced (strongly, 
in the case of the outermost layer) with respect to the bulk, for both types of systems, despite its 
reduced dimensionality in metals. 

PACS numbers: 61.25.Mv, 64.70.Fx, 71.15.Pd 
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I. INTRODUCTION. 

Surface layering at metallic liquid-vapor interfaces 
was suggested in the early eighties by Monte Carlo 
simulations of alkali metal liquid surfaces 0. Subse- 
quently, more elaborate Monte Carlo techniques (see 
[3| and references therein), experimental measurements 
13, 0, 0, 0, 0, IS S ^1 and recent ab initio molecular dy- 
namics (MD) simulations [H El E El El , have estab- 
lished that the ionic density profile of liquid metals across 
the interface shows oscillations that decay into the bulk 
liquid after, on average, three to four layers. The use of 
ab initio MD simulations in these studies is important for 
two reasons. First, being ab initio, the valence electrons 
and the ions are treated on the same footing, with the 
ions reacting consistently to the large spatial variations 
undergone by the electronic density when moving from 
the liquid to the vapor. Second, being MD simulations, 
a study of the motion of the atoms can be carried out di- 
rectly and changes across the interface can be analyzed. 
This type of study has not yet been undertaken, as far 
as we know, basically because most of the simulations of 
liquid metallic surfaces were carried out using the Monte 
Carlo method, which gives no information about the dy- 
namic properties of the system studied. The main aim 
of this paper is to fill this gap, reporting a study of the 
atomic motion in the liquid-vapor interfaces of Li, Na, K, 
Rb, Cs, Mg, Ba, Al, Tl, Si, and the NasKy alloy, using 
slabs with 2000 atoms. 

Layering also appears for metals at the solid-liquid 
interface [iGj, for liquids in contact with a hard wall 
[TtL IT^ . or in confined geometries 0|. Moreover, for 
ultrathin films of large molecules in contact with a Si 
wall ^0] and for confined dusty plasma liquids E3 it has 
been possible to follow experimentally the molecular mo- 
tion, whereas in the solid-liquid interface ab initio MD 



simulations allowed to study the atomic diffusion at the 
interface (16] . In all cases reduced diffusion with respect 
to the bulk was found. 

These oscillating profiles are in contrast with those of 
the liquid-vapor interfaces of one-component dielectrics, 
like water or Lennard-Jones (LJ) systems, or those of 
liquid-liquid interfaces of immiscible mixtures, which 
show monotonic profiles with no layering. Taking the 
archetypical example of the one-component LJ system, 
we have found in the literature a large amount of com- 
puter simulation studies of the coexisting densities, the 
interfacial width and the surface tension, as a function 
of the temperature, the particular LJ model used (trun- 
cated, truncated and shifted, full potential) and the lat- 
eral area simulated (see, for instance, and references 
therein). However, to our knowledge, an analysis of the 
atomic motion in the interface has not been performed 
yet. Surprisingly, the situation is somewhat different for 
more complex systems, like water, or mixtures, where 
some (scarce) studies of atomic diffusion in the interfa- 
cial region have indeed been performed |2^, "2?, "251 . 
In the case of mixtures, some anisotropy in diffusion was 
detected, either due to an increase in the diffusion coef- 
ficient parallel to the interface, Dt, with respect to the 
bulk value 22, 23], or to a decrease of the diffusion co- 
efficient normal to the interface, D^, Following a 
thorough analysis of the effect of the inhomogeneity of 
the interface on the values (and even the definition) of 
-Dat, molecular dynamics simulations of the liquid- vapor 
interface of water |2l revealed an increase in Dt (3.5- 
fold) and also in Dn (2-fold), which was attributed to the 
reduction of the number of hydrogen bonds in the inter- 
face. Moreover, it was argued that the same qualitative 
behavior would occur for other simple liquids 2^^ . 

Although surface layering is common to the liquid- 
vapor interface of metals and the systems mentioned 
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above, there are also important differences. Even though 
the atoms in the outer layer of the liquid- vapor interface 
rarely leave the liquid surface (especially at temperatures 
near their triple points), they are not geometrically con- 
fined. Moreover, in liquids confined, or on a wall, or at 
the interface with their solid phase, a strong interaction 
between the substrate and the first layers of the liquid 
has an important influence on the structure and dynam- 
ics of the liquid, leading to strong oscillations in the pro- 
file and to reduced diffusion. It is also worth analyzing 
if the metallic liquid-vapor interface, in particular the 
outer layer, behaves like a quasi two-dimensional system. 
There are some liquid alloys (Ga with small amounts of 
Tl, Pb or Bi) where the minority component, which has 
a high melting temperature and segregates to the sur- 
face . di splays this kind of behavior (see, e. g., references 
|2a. I27l l2 q | ) . For one component metals or other type of 
alloys, a possible measure of the two-dimensional char- 
acter of the outer layer could be some distinctive feature 
of the density distribution function (DDF) of the time of 
residence of the atoms in that layer, for instance its mean 
value. 



II. SIMULATION DETAILS. 

Because of the the periodic boundary conditions used 
in most simulation methods, a slab geometry is usually 
adopted in studies of the liquid- vapor interface, with two 
free surfaces perpendicular to one of the axes, taken here 
as the z axis. These two surfaces should be well separated 
to minimize interaction between them. For the metallic 
systems we have performed ab initio molecular dynamics 
simulations, whereas for the LJ system standard classical 
molecular dynamics have been used. 

Ab initio simulations based on the Kohn-Sham formu- 
lation of Density Functional Theory (DFT) j pose 
huge computational demands, and the only metallic liq- 
uid surfaces studied so far. Si ^llj and Na jl^j, used small 
samples (96 and 160 particles respectively) where the two 
surfaces are rather close (16 and 20 A respectively). Or- 
bital free ab initio molecular dynamics (OFAIMD) simu- 
lations reduce somewhat these computational demands, 
by returning to the original Hohenberg-Kohn formula- 
tion of DFT [s^l, and adopting an explicit but approx- 
imate density functional for the electron kinetic energy 
so that the whole energy is a functional of the electron 
density. OFAIMD simulations of li quid metal surfaces 
have recently been performed [T^l 0, 15] for Li, Na, Mg, 
Al, Si and Na3K7 using 2000 particles and for Li4Na6 
using 3000 particles, which led to simulation boxes big 
enough to reliably represent a macroscopic interface. De- 
tails about the formalism and the electron-ion interac- 
tions can be found elsewhere [s^ • Even though the 
OFAIMD method has been applied successfully to bulk 
metals and alloys, it might be argued that this is not a 
sufficient validation of the method, as surfaces are dif- 
ferent from bulk systems. However, it must be stressed 



that the OFAIMD studies of Na and Si, produced re- 
sults very similar to those obtained by the (in princi- 
ple) more accurate Kohn-Sham ah initio simulations: the 
wavelength of the oscillations in the profiles, which is re- 
covered exactly, the number of nearest neighbors of a Si 
particle across the interface, and the surface tension of 
liquid Na, are all well reproduced by the OFAIMD ap- 
proach. Further confidence in the ability of the methd to 
tackle metallic surfaces can be obtained from studies for 
finite systems, where the surface plays an essential role. 
For instance, a long standing and previously unexplained 
anomalous variation of the melting temperatures of Na 
clusters with size |33| has been reproduced and rational- 
ized for the first time in terms of the surface geometry 
and stability using the same OFAIMD method as in 
this paper. Even furher confidence in the capabilities of 
the method, concerning specifically semiintinite surfaces, 
can be gained from preliminary results ^35] obtained for 
the temperature dependent surface relaxation of Al(llO) 
and Mg(lOTO), which reproduce qualitatively both exper- 
imental data and previous Kohn-Sham calculations. 

The LJ system has been simulated at a reduced tem- 
perature T* = ksT/e = 0.73, which is near its triple 
point, in order to compare with the metallic systems in 
a similar thermodynamic state. We have made simula- 
tions with 1960 particles, similar to the metals, and also 
with 15680 particles in order to reassess the results for 
the smaller system. The lateral side for the small system 
is 11.75cr, which is already large enough to suppress the 
unrealistic oscillatory behaviour of the interfacial prop- 
erties due to periodic boundary conditions [s^. For the 
large system the lateral side is doubled (23.5(t) so the ex- 
pected effects on the interfacial properties are only those 
of the enhanced capillary waves which should produce a 
wider interface than in the small system. 

For all systems, the number of NVE equilibrium con- 
figurations used for averaging ranged between 20000 and 
30000. 



III. LAYERS DEFINITION. 

Figure ^ shows the ionic density profile obtained for 
liquid Si and the partial and total ionic density profiles 
for the liquid alloy NasKr, while that of the LJ system is 
depicted in figure [3 The layers where the ionic motion 
was studied, which for metals were located between con- 
secutive minima in the total ionic density profile, are also 
indicated. The outermost layer (numbered 1) comprises 
from the last minimum to the infiection point of the de- 
caying profile. For layers 1, 2 and 3 the results of the 
two layers on opposite sides of the slab were averaged. 
Further structural details of the interfaces of these sys- 
tems will be given elsewhere [s^, but we note here that 
the relative amplitudes of the outermost oscillation for 
the alkalis and alkaline-earths are rather similar (max- 
ima from 1.00 to 1.13, minima from 0.80 to 0.88) while 
they increase significantly for the systems of valence 3 
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FIG. 1: Ionic density profiles normal to the liquid-vapor in- 
terfaces of Si and NasKy. The different regions where the 
atomic motion has been studied are also shown. 



FIG. 3: Time evolution of the z coordinate of an Al atom 
and its projected trajectory onto the xz (full line) and the yz 
(dotted line) planes. 



and 4 (maxima 1.29, 1.47 and 1.56, minima 0.73, 0.69 
and 0.54 for Al, Si and Tl respectively). 

The width of the different layers is very small, of the 
order of one atomic diameter. Therefore it is not appro- 
priate to talk about diffusion along the direction normal 
to the interface within a layer, as there is no room for 
a particle within the layer to reach that type of motion. 
Figure Ol shows the time evolution of the z coordinate of 
one particular Al atom and an oscillatory-type motion 
within the layers can be seen followed by jumps into ad- 
jacent layers. The trajectory of the particle, projected 
onto the xz and yz planes is also shown in the figure. 
To provide a meaningful comparison among the results 
for different systems, with diverse masses and potentials, 
and consequently different characteristic times, a refer- 
ence time, r, has been defined for each system which is 
to be used as a time unit (see table . Specifically, r 
has been taken as the first value of t for which the mean 
square displacement (MSD) in the bulk liquid has an in- 
flection point, which is an indication of the atomic motion 
changing from free-particle to diffusive like behavior. 
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The ionic density profile for the large LJ system is 
shown in figure HJ For this type of system, the definition 
of "layers" in order to analyze atomic motion is basically 
arbitrary. We have chosen to define them in a similar 
way as we have done for metals in order to make the 
comparison as fair as possible. First we have fitted the 
simulation results to an error function profile, 

P(^) ^\{Pt + Pv) + ^{pe- p^,)erf(^/¥(zo - z)/w) 

(where pi and py are the coexisting densities of the liquid 
and the vapor, zq is the position of the interface, which 
coincides with the inflection point of the profile, and w 
its width), which appears to be more adequate than the 
usual hyperbolic tangent one [sH. The widths for the 
small and large LJ systems are w = 2.25a and w = 2 Ala 
respectively, the latter being larger as expected from the 
increased number of capillary waves. In the following 
we will present the results obtained for the large system, 
which, appart from this interfacial width, are practically 
coincident with those of the small system. Layer number 
1 is then defined as a slice of width w /2 from the inflection 
point towards the liquid. Further slices of width 1x1/2 
towards the liquid are then taken as layers 2 and 3. Note 
that, defined this way, the layers of the LJ system are also 
very thin (1.235(t), so we are in a similar situation as in 
metals, where we consider that diffusion is only possible 
within the layers parallel to the interface. The reference 
time r for the LJ system is finally defined in the same 
way as for metals, and is also shown in table 



FIG. 2: Ionic density profile normal to the liquid- vapor inter- 
faces of the LJ system at T* = 0.73, and the "layers" defined 
for this interface. Only half the slab is shown. 



IV. RESULTS AND DISCUSSION. 

In order to analyze the diffusive motion within the lay- 
ers, the MSD along the x and y directions has been cal- 
culated for those particles that remain inside the layer 
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FIG. 4: Distribution function of the residence time, p(t), in 
layer 1 (full line) and layer 3 (dotted line) for liquid Ba. The 
inset shows the integrated probability distribution function, 
P{t), together with the 20 % level used to define t2o- The 
dashed line denotes layer 2. 

for a large enough time (tm), a criterion somewhat dif- 
ferent from, althou gh i n the same spirit as others used in 
previous works p3.l25| . The selection of tm is somewhat 
delicate, as tm must be large enough so that diffusive 
motion has already set in, but the number of particles 
that remain in a layer for at least tm decreases as tm in- 
creases, leading to poorer statistics. A good compromise 
was found to be tm — ^20 j defined so that the probability 
of residence in the layer for times larger than t2o is 20 %. 
All the data shown below were obtained with this cri- 
terion, but checks showed that other reasonable choices 
for tm led to statistically equivalent results. In order to 
compute ^20, and also to characterize the motion along 
the z-direction, the residence time of particles in each of 
the layers was studied. 



A. Motion perpendicular to the interface. 

For all the systems considered, including the metals 
and the LJ system, the DDF of the residence time, p{t), 
has been computed for the three outer layers (the case of 
Ba is shown in figure 0}, and some common features can 
be identified. The maximum of p{t) occurs at very small 
i- values, (« r), suggesting that most of the particles at- 
tempting to enter the layer are bounced back. Also, the 
DDF exhibits a long time tail, which gives rise to mean 
residence times (MRT), iav, which increase when mov- 
ing from the bulk liquid towards the outer layers (see 
table nj. In the LJ system this increase is rather small, 
amounting to approximately a 10 %. For the metals, 
however, the MRT at the outermost layer is strongly en- 
hanced, with increases between 75 % (for Na) and 246 % 
(for Tl), while the values taken at layer 2 show a more 
modest increase of around 10 %. This large discrepancy 



FIG. 5: Mean square displacements for layers 1 and 2 and the 
slab center for liquid K. 



between the behavior of the LJ system and that of met- 
als precludes any simple geometrical interpretation of the 
results and underlines that the change of interactions at 
the liquid- vapor interface of metals does have an impor- 
tant influence on some dynamic properties of the surface. 
Moreover, these results indicate that the outermost layer 
in metals behaves more bidimensional-like than the inner 
ones. It is interesting to note that the MRT is around 20r 
for the alkalis and alkaline-earths, whereas it increases in 
Al (w 40r) and Tl (« 60t), but not in Si. 

In order to rationalize these differences among metals, 
we consider how the atoms move from one layer to an- 
other. The minima in the density profiles suggest the 
existence of an interlayer potential barrier, which would 
be higher the lower the minimum of the density profile. 
Therefore this barrier would increase significantly from 
the alkalis and alkaline-earths to Al, to Si, and finally to 
Tl. The probability of overcoming this barrier would be 
related both to its height and to the frequency with which 
the atoms attempt to cross. This frequency has been es- 
timated as proportional to the Einstein frequency along 
the z direction, fi^, which is shown in table together 
with the corresponding frequency in the direction par- 
allel to the interface, fix, for completeness. These have 
been obtained from a short time (up to t = t) expansion 
of the corresponding MSD 39]. Wc find, in terms of r, 
a rather universal value for all the systems except for Si, 
whose more open structure leads to higher Einstein fre- 
quencies. Therefore, the increase in the MRT of Al and 
Tl is related to the increased barrier height, whereas in 
Si the increased barrier height is counterbalanced by an 
increased Einstein frequency, leading to MRT similar to 
those of alkalis and alkaline-earths. 

The absence of well defined layers in the LJ profile 
suggests that the motion perpendicular to the interface 
is much easier than in metals, as no barriers are present, 
and therefore the position of our arbitrarily defined "lay- 
ers" has very little influence on the MRT. 
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System 


r 


Layer 1 


Layer 2 


Layer 3 


Layer 1 


Layer 2 


Layer 3 


L 

Layer 1 


>T / -Dcent 

Layer 2 


2r 

Layer 3 


S 

Layer 1 


Layer 2 


Center 


LJ 


0.167 


11.1 


10.6 


10.5 


18.7 


17.7 


17.2 


1.94 


1.13 


0.96 


-(-) 


-(-) 


— 


Li 


0.038 


24.2 


14.2 


13.2 


32.4 


19.7 


17.4 


2.01 


1.19 


1.09 


1.27(1.18) 


1.58(1.51) 


1.56 


Na 


0.105 


17.6 


11.0 


10.1 


26.3 


16.6 


14.6 


2.16 


1.18 


1.12 


1.37(1.29) 


1.65(1.59) 


1.60 


K 


0.170 


19.5 


11.6 


10.8 


26.8 


16.8 


15.0 


2.19 


1.21 


1.08 


1.35(1.28) 


1.64(1.57) 


1.59 


Rb 


0.284 


19.5 


11.3 


10.6 


27.7 


16.9 


15.6 


2.10 


1.22 


1.08 


1.36(1.29) 


1.67(1.59) 


1.61 


Cs 


0.352 


20.3 


12.0 


10.6 


32.2 


17.0 


15.3 


1.81 


1.13 


1.07 


1.41(1.35) 


1.73(1.65) 


1.67 


Mg 


0.056 


22.5 


12.0 


11.1 


33.7 


17.7 


15.2 


2.27 


1.29 


1.09 


1.33(1.33) 


1.65(1.58) 


1.61 


Ba 


0.158 


24.6 


13.5 


12.8 


25.1 


16.5 


14.5 


2.02 


1.22 


1.08 


1.30(1.28) 


1.63(1.56) 


1.55 


Al 


0.045 


37.3 


17.6 


14.4 


68.0 


25.6 


19.3 


2.04 


1.40 


1.33 


1.32(1.46) 


1.66(1.62) 


1.65 


Tl 


0.187 


64.7 


26.1 


18.7 


108.2 


39.4 


27.0 


1.24 


1.00 


1.00 


1.32(1.46) 


1.65(1.60) 


1.62 


Si 


0.058 


21.4 


9.5 


7.8 


37.4 


16.3 


13.3 


1.24 


1.00 


1.00 


1.68(1.87) 


2.05(1.94) 


1.97 


NaQNagKr 


0.126 


18.7 


10.5 


17.4 


14.7 


1.28 


1.05 


1.61(1.57) 


1.61 


KONagKr 


0.176 


18.9 


10.2 


9.8 


29.4 


15.2 


14.4 


2.24 


1.25 


1.06 


1.39(1.28) 


1.71(1.64) 


1.67 



TABLE I: Reference time, r, mean residence time, fav, twenty percent time (see text), t2o, ratio between Dt and the diffusion 
coefficient in the center of the slab, Dccntcr, and Einstein frequencies for the systems considered and the different regions. The 
units of r for the LJ system are standard reduced units, and picoseconds for all the other systems. 



Similar trends are also obtained for the liquid NasKy 
alloy, although surface segregation leads to an outer layer 
of almost pure K so that Na-related quantities are really 
irrelevant. Indeed the outermost Na layer comprises re- 
gions 1 and 2, and therefore the numbers shown in table 
IHfor Na in the alloy span columns 1 and 2 together. 

I 



B. Motion parallel to the interface. 

The ^20 times have been used to analyse the atomic 
motion within the layers. The MSD have been evaluated 
up to 0.8 X t20; in order to allow an adequate averaging 
over the time origins. Figure|Sldepicts the MSD for layers 
1 and 2 and for the slab center in liquid K. There is a 
clear increase in the slope as the interface is approached. 
This is common to all the systems, including the LJ one, 
and is quantified in table^which shows the ratio between 
the diffusion coefficient parallel to the interface, Dt, in 
the different layers and the bulk diffusion coefficient in 
the slab center (where diffusion is isotropic). 

We attribute this ~ 100 % increase in the diffusion at 
the outermost layer to the reduced coordination of the 
atoms in the interface. In all the pure metals, except Si, 
the number of nearest neighbors is w 12 at the center of 
the slab and reduces to ~ 8 at the outer layer. In the 
LJ system the cordination number again decreases from 
« 12 to ~ 9. For Si the coordination number is « 6 at 
the center and « 4.5 at the outer layer. This smaller 
decrease is reflected in a smaller enhancement of the dif- 
fusion coefficient at the interface of only 24 %. A similar 
explanation also holds for the NasKy alloy. The total 
number of neighbors around a K (Na) atom is roughly 



The integral P{t) = p{u)du, gives the probability 
that having entered the layer, a particle remains there 
longer than t, and is used to obtain i20 (see the inset of 
figure Values for ^20 shown in table correlate rather 
well with the MRT. 



I 

constant at « 12.5 (10.0) up to regions 1, where it sharply 
decreases. Consequently, Dt for K increases rapidly in 
the outermost layer, whereas the value for Na in its outer- 
most layer (regions 1 and 2) does not change significantly 
because very few Na atoms in regions 1 are affected by 
this decrease of coordination. 

The case of Tl deserves special attention because the 
number of neighbors decreases from 12 to 8 but Dt in- 
creases by only 24 % in the outermost layer. We attribute 
this modest increase in Dt in this strongly layered sys- 
tem to a large influence of the second layer on the atoms 
of the outermost one, an effect similar to that exerted 
by a wall on a liquid which leads to strong layering and 
reduced diffusion. 



V. CONCLUSSIONS. 

In summary, we have analyzed through ab initio simu- 
lations the atomic motion in the liquid- vapor interfaces of 
several metals and compared it with that of a L J system 
in a similar thermodynamic state. Although the layered 
structure is similar to other systems such as liquids on 
a solid or in confined geometries, the dynamic behavior 
within the layers is much more similar to the liquid- vapor 
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interface of dielectrics, like the LJ system or water, show- 
ing an enhanced diffusion in the parallel direction at the 
interface, which is attributed to the reduced coordina- 
tion of an atom which favours the transverse movement. 
In the perpendicular direction, the layers are too thin to 
regard the motion as atomic diffusion, and instead the 
MRT associated to each layer has been determined. The 
dynamic behavior of metallic systems along this direc- 
tion is far different from that of LJ systems. The value 



of the MRT for metals is clearly larger at the outermost 
layer, contrary to the case of the LJ system, and increases 
significantly for the polyvalent metals with closed struc- 
tures. Finally, the reference time r is found to be an 
excellent time unit, since it reveals universal values in 
several dynamic properties of different systems. 

We acknowledge the financial support of the DGICYT 
(MAT2005-03415) and the EU FEDER program. 
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